%% process_pde_stokes_s
function pde = process_pde_stokes_s(index)
switch index
    case 1
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y) x.^2 .* y.^2 + exp(-y);
        pde.u1dx = @(x,y) 2*x.*y.^2;
        pde.u1dy = @(x,y) 2*y.*x.^2 - exp(-y);
        pde.u2 = @(x,y) -2/3*x.*y.^3 + 2 - pi*sin(pi*x);
        pde.u2dx = @(x,y) - pi^2*cos(pi*x) - (2*y.^3)/3;
        pde.u2dy = @(x,y) -2*x.*y.^2;
        pde.p = @(x,y) -(2-pi*sin(pi*x)).*cos(2*pi*y);
        pde.pdx = @(x,y) pi^2*cos(pi*x).*cos(2*pi*y);
        pde.pdy = @(x,y) -2*pi*sin(2*pi*y).*(pi*sin(pi*x) - 2);
        pde.f1 = @(x,y) -2*nu*x.^2 - 2*nu*y.^2 - nu*exp(-y) + pi^2*cos(pi*x).*cos(2*pi*y);
        pde.f2 = @(x,y) 4*nu*x.*y - nu*pi^3*sin(pi*x) + 2*pi*(2-pi*sin(pi*x)).*sin(2*pi*y);
        pde.domain = [0, 1, -0.25, 0];
        pde.bdary = [0, 0, 0, 0];
    case 2
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y) pi*sin(2*pi*y).*sin(pi*x).^2;
        pde.u1dx = @(x,y) pi.^2.*cos(x.*pi).*sin(x.*pi).*sin(y.*pi.*2.0).*2.0;
        pde.u1dy = @(x,y) pi.^2.*cos(y.*pi.*2.0).*sin(x.*pi).^2.*2.0;
        pde.u2 = @(x,y) -pi*sin(2*pi*x).*sin(pi*y).^2;
        pde.u2dx = @(x,y) pi.^2.*cos(x.*pi.*2.0).*sin(y.*pi).^2.*-2.0;
        pde.u2dy = @(x,y) pi.^2.*cos(y.*pi).*sin(x.*pi.*2.0).*sin(y.*pi).*-2.0;
        pde.p = @(x,y) 10*cos(pi*x).*sin(pi*y);
        pde.pdx = @(x,y) pi.*sin(x.*pi).*sin(y.*pi).*-1.0e+1;
        pde.pdy = @(x,y) pi.*cos(x.*pi).*cos(y.*pi).*1.0e+1;
        pde.f1 = @(x,y) -nu.*(pi.^3.*cos(x.*pi).^2.*sin(y.*pi.*2.0).*2.0-pi.^3.*sin(x.*pi).^2.*sin(y.*pi.*2.0).*6.0)-pi.*sin(x.*pi).*sin(y.*pi).*1.0e+1;
        pde.f2 = @(x,y) pi.*cos(x.*pi).*(cos(y.*pi).*5.0-nu.*pi.^2.*sin(x.*pi).*6.0+nu.*pi.^2.*cos(y.*pi).^2.*sin(x.*pi).*8.0).*2.0;
        pde.domain = [0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 3
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y) 10.*x.^2.*y.*(2.*y - 1).*(x - 1).^2.*(y - 1);
        pde.u1dx = @(x,y) x.*y.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*2.0e+1;
        pde.u1dy = @(x,y) x.^2.*(x-1.0).^2.*(y.*-6.0+y.^2.*6.0+1.0).*1.0e+1;
        pde.u2 = @(x,y) -10.*x.*y.^2.*(2.*x - 1).*(x - 1).*(y - 1).^2;
        pde.u2dx = @(x,y) y.^2.*(y-1.0).^2.*(x.*-6.0+x.^2.*6.0+1.0).*-1.0e+1;
        pde.u2dy = @(x,y) x.*y.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*-2.0e+1;
        pde.p = @(x,y) (20.*x - 10).*(2.*y - 1);
        pde.pdx = @(x,y) y.*4.0e+1-2.0e+1;
        pde.pdy = @(x,y) x.*4.0e+1-2.0e+1;
        pde.f1 = @(x,y) (y.*2.0-1.0).*(nu.*y-nu.*x.^2.*3.0+nu.*x.^3.*6.0-nu.*x.^4.*3.0-nu.*y.^2+nu.*x.*y.^2.*6.0+nu.*x.^2.*y.*6.0-nu.*x.^2.*y.^2.*6.0-nu.*x.*y.*6.0+1.0).*2.0e+1;
        pde.f2 = @(x,y) (x.*2.0-1.0).*(-nu.*x+nu.*x.^2+nu.*y.^2.*3.0-nu.*y.^3.*6.0+nu.*y.^4.*3.0-nu.*x.*y.^2.*6.0-nu.*x.^2.*y.*6.0+nu.*x.^2.*y.^2.*6.0+nu.*x.*y.*6.0+1.0).*2.0e+1;
        pde.domain = [0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 4
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y) 5*(x.^2 - 2*x.^3 + x.^4).*(2*y - 6*y.^2 + 4*y.^3);
        pde.u1dx = @(x,y) x.*y.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*2.0e+1;
        pde.u1dy = @(x,y) x.^2.*(x-1.0).^2.*(y.*-1.2e+1+y.^2.*1.2e+1+2.0).*5.0;
        pde.u2 = @(x,y) -5*(y.^2 - 2*y.^3 + y.^4).*(2*x - 6*x.^2 + 4*x.^3);
        pde.u2dx = @(x,y) y.^2.*(y-1.0).^2.*(x.*-1.2e+1+x.^2.*1.2e+1+2.0).*-5.0;
        pde.u2dy = @(x,y) x.*y.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*-2.0e+1;
        pde.p = @(x,y) 5*(2*x - 1).*(2*y - 1);
        pde.pdx = @(x,y) y.*2.0e+1-1.0e+1;
        pde.pdy = @(x,y) x.*2.0e+1-1.0e+1;
        pde.f1 = @(x,y) y.*2.0e+1-nu.*(x.^2.*(y.*2.4e+1-1.2e+1).*(x-1.0).^2.*5.0+y.*(x.*-6.0e+1+x.^2.*6.0e+1+1.0e+1).*(y.*-3.0+y.^2.*2.0+1.0).*2.0)-1.0e+1;
        pde.f2 = @(x,y) x.*2.0e+1+nu.*(y.^2.*(x.*2.4e+1-1.2e+1).*(y-1.0).^2.*5.0+x.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-6.0e+1+y.^2.*6.0e+1+1.0e+1).*2.0)-1.0e+1;
        pde.domain = [0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    otherwise
        error("Invalid pde index.");
end
fprintf("pde index: %d\n", pde.index);
fprintf("pde domain: xmin,xmax,ymin,ymax\n");
disp(pde.domain);
fprintf("pde boundary: bottom,right,top,left (0=Dirichlet,1=Neumann)\n");
disp(pde.bdary);
if isfield(pde, "bdary_p")
    fprintf("pde boundary of p: (0=Dirichlet,1=Neumann)\n");
    disp(pde.bdary_p);
end
end
